function log_like = LLcorr_probit(param)

%% GLOBAL

global Y T1 T2 SM SR Draws M Xc Gmax Gmin or orm ll


%% PARAM def 


cm = param(1);
cr = param(2);
phim = param(3);
phir = param(4);
sigem = param(5);
siger = param(6);
cgpa = param(7);
phigpa = param(8);
sigegpa = param(9);

c1 = param(10);
c2 = param(11);
b1 = param(12);
b2 = param(13);
phi1 = param(14);
phi2 = param(15);
gamma1 = param(16);
gamma2 = param(17);
sigtheta = param(18);
sige12 = param(19);
cy = param(20);
by = param(21);

sige1 = 1;
sige2 = 1;
sigy = 1;

c1d = param(22);
c2d = param(23);
phi1d = param(24);
phi2d = param(25);
b1d = param(26);
b2d = param(27);





%%
sigtheta = abs(sigtheta);
theta = Draws*sigtheta;

Lsm = normalpdf(repmat(SM,1,M)-cm-phim*theta,0,sigem);
Lsr = normalpdf(repmat(SR,1,M)-cr-phir*theta,0,siger);



if abs(sige12)>1
    sige12=sige12/(abs(sige12)+0.0001);
end;

Tt1=repmat(T1,1,M); Tt2=repmat(T2,1,M);

w=[reshape((2*Tt1-1).*(c1+phi1*theta+b1*repmat(Xc,1,M)+c1d*or+phi1d*or.*theta+b1d*or.*repmat(Xc,1,M)),[],1),reshape((2*Tt2-1).*(c2+c2d*orm+phi2*theta+phi2d*orm.*theta+b2*repmat(Xc,1,M)+b2d*orm.*repmat(Xc,1,M)),[],1)];
t1=reshape(Tt1,[],1); t2=reshape(Tt2,[],1);
Prob=((t1.*t2).*mvncdf(w,[0,0],[1,sige12;sige12,1]))...
    +((t1.*(1-t2)).*mvncdf(w,[0,0],[1,-sige12;-sige12,1]))...
    +(((1-t1).*t2).*mvncdf(w,[0,0],[1,-sige12;-sige12,1]))...
    +(((1-t1).*(1-t2)).*mvncdf(w,[0,0],[1,sige12;sige12,1]));

Lt12=reshape(Prob,[],M);
Lgpa = normcdf(cgpa+phigpa*theta-Gmax,0,sigegpa).*(repmat(Xc,1,M)>=Gmax)+(repmat(Xc,1,M)<Gmax & repmat(Xc,1,M)>Gmin).*normalpdf(repmat(Xc,1,M)-cgpa-phigpa*theta,0,sigegpa)+(repmat(Xc,1,M)==Gmin).*(1-normcdf(cgpa+phigpa*theta-Gmin,0,sigegpa));

%% Bias

yind = normcdf(cy+theta+by*repmat(Xc,1,M));
yind = yind+0;
bias1 = repmat(T1,1,M)-yind;
bias2 = repmat(T2,1,M)-yind;


Ly = (normcdf(cy+theta+gamma1*bias1+gamma2*bias2+by*repmat(Xc,1,M), 0, sigy).^repmat(Y,1,M)).*((1-normcdf(cy+theta+gamma1*bias1+gamma2*bias2+by*repmat(Xc,1,M), 0, sigy)).^(1-repmat(Y,1,M)));

%% Individual Likelihood Contribution
Lm=(((Lsm.*Lsr).*Lt12).*Ly).*Lgpa;


%% log likelihood
ll = log(mean(Lm,2));
if abs(param(18))<0.1;
    log_like=1e5;
else
log_like = -sum(log(mean(Lm,2)));
end;

